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ABSTRACT 

Aims. We present a new formal solution of the Lagrangian equation of radiative transfer that is useful in solving the equation of 
radiative transfer in the presence of arbitrary velocity fields. 

Methods. Normally a term due to the inclusion of the wavelength derivative in the Lagrangian equation of radiative transfer is 
associated with a generalised opacity. In non-monotonic velocity fields, this generalised opacity may become negative. To ensure that 
the opacity remains positive, this term of the derivative is included in the formal solution of the radiative transfer problem. 
Results. The new definition of the generalised opacity allows for a new solution of the equation of radiative transfer in the presence of 
velocity fields. It is especially useful for arbitrary velocity fields, where it effectively prevents the occurrences of negative generalised 
opacities and still allows the explicit construction of the A-operator of the system needed for an accelerated A-iteration. We performed 
test calculations, where the results of old, established solutions were compared with the new solution. The relative deviations never 
exceeded 1 % and so the new solution is indeed suitable for use in radiative-transfer modelling. 

Non-monotonic velocity fields along photon paths frequently occur in three-dimensional hydrodynamical models of astrophysical 
atmospheres. Therefore, the formal solution will be of use for multidimensional radiative transfer and has immediate applications in 
the modelling of pulsating stars and astrophysical shock fronts. 
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1. Introduction 

The Lagrangian equation of radiative transfer for moving atmo- 
spheres (Mihalas 1980) has been successfully solved for differ- 
ent astrophysical atmospheres, in particular supernova and nova 
atmospheres, as well as stellar atmospheres with winds. 

The inclusion of velocity fields in the equation of transfer 
gives rise to two terms, one which is proportional to the spe- 
' cific intensity. This term is of the same form as the term aris- 
' ing from the opacity. Thus, the term from the velocity field can 
be included in the definition of a generalised opacity (see for 
instance Hauschildt 1992). For positive, velocity gradients, the 
contribution of the velocity field to the generalised opacity is 
always positive. In the case of negative, velocity gradients, the 
velocity field contribution may become dominant and produce a 
negative, generalised opacity. Thus, the method of solution must 
avoid exponentially increasing terms in the formal solution be- 
cause such terms make the numerical computation unstable. 

We show how to avoid this problem by not introducing a 
generalised opacity but instead including the velocity field in 
the formal solution of the radiative transfer. Negative opacity 
will occur in non-monotonic velocity fields, thus any numeri- 
cal scheme for the formal solution must be able to handle both 
positive and negative velocity gradients. Numerically, the for- 
mal solution of the radiative transfer problem is then no longer 
an initial-value problem but a boundary-val ue problem. A formal 
solution of this situation was described by iBaron & Hauschildl 



(12004 ') and the revised formal solution in this paper is derived 
from that solution. In addition, the formal solution incorporates 
the two discretisation schemes fo r the waveleng th derivative and 
their mixing described by H auschildt & BaronI (,2004). It it im- 
portant to note that the new formal solution must be consistent 
with the solution to the scattering problem, e.g. an Accelerated 
A-iteration (ALI) (.Scharmer 1981) as used in this work. 



An immediate application of this new solution is to both 
one-dimensional model atmospheres of pulsating giant stars and 
shock fronts of all kinds, such as the interaction of a supernova 
envelope with the interstellar medium producing both forward 
and reverse shocks. In the course of multidimensional modelling 
of astrophysical flows, the use of these solutions will become un- 
avoidable since arbitrary velocity fields occur frequently, due to 
e.g. Rayleigh-Taylor instabilities. In fact it is impossible to per- 
form comoving radiative transfer in arbitrary velocity fields by 
means of a characteristic method without the use of a formal so- 
lution that totally avoids the occurrence of negative generalised 
opacities. 



We derive the revised formal solution with a different defini- 
tion of the opacity in Sect.|2l In Sect.[3j we demonstrate that this 
solution produces results consistent with established solutions 
and present brief conclusions in Sect.|4] 
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2. The construction of the formal solution 

We solve the equation of radiative transfer by following the so- 
lution along characteristic rays. In this form, the transport of the 
specific intensity / is described along a photon path (with path 
length s) and is independent of the dimensionality of the under- 
lying atmosphere: 



d^ 
di 



■ aA- 



dA ' 



where is the emissivity and x is the total extinction opacity. 
The dependence on the comoving wavelength has been explic- 
itly separated from the differential A giving rise to the 
terrr(3 that can be incorporated into a generalised opacity. The 
term = |^ is due to the wavelength coupling and is, for in- 
stance, given in spherical symmetry by: 



\ or r 



(2) 



where (3 - v/c, y - Ij yj\ — 0^, and r and yu - cos j? are 
the phase-space coordinates (wavelength is of course a phase- 
space coordinate). We note that it is the sign of the wavelength- 
coupling term, and not the sign of the velocity gradient, that de- 
termines the upwind sense of the wavelength derivative. 

The form of Eq.[T]is the same for all geometries and dimen- 
sionalities. However, the form of the differential 4- and the cou- 

' as 

pling term ax will differ 

To formally solve Eq.[Tl the wavelength derivative must 
be dis cretised. We adopt the notation of iHauschildt & BaronI 
(|2004|) . where two possible discretisations of the wavelength 
derivative are combined in a Crank-Nicholson-like scheme. 

Since we allow for arbitrary velocity fields, the formal so- 
lution must allow for an arbitrary sense of wavelength deriva- 
tive that is, the upwind direction may correspond either to longer 
or shorter wavelengths. To insure the stability of the discre tisa- 
tion, local upwind schemes are introduced dBaron & Hauschildtl 
l2004i) depending on the coupling term a^- We assume a sorted 
wavelength grid < Ai < Ai+\, and the wavelength depen- 
dence is represented by the wavelength index. The wavelength 
derivative can then be written as: 



dA 
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where the p' are defined as: 
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It should be noted that the p' depend not only on wavelength but 
also on spatial position since the sign of the coupling term 
may change along the characteristic. 

For mixing parameters of ^ e [0, 1] for the two different 
discretisations, the equation of radiative transfer then reads: 



As 



Ii-a,{4 + ^ p|) Ii-^ai [pi Ii^i + pIIm) 



-[l-^]ai [p, h-i + p\li + pIIm) 



(4) 



To solve Eq.lHfor monotonic velocity fields, it is customary 
to define a generahsed opacity: 



Xi^Xi+ai{4 + ^ p^^ 



(5) 



This approach in general, however, often fails because the 
generalised opacity in Eq. |5]may become negative if: 



(1) - 4ai > xi + ai^ p\ 



(6) 



(Note that the coefficient always has the same sign as ai). For 
strong wavelength couplings — as for instance in shock flows 
— the condition |6] is easily fulfilled. To avoid negative optical 
depths along the characteristics, the generalised opacity must be 
defined differently. In the ^ - I case, negative opacities could 
be eliminated by adopting a fine wavelength sampling anywhere 
apart from at the boundaries where the vanish. In principle, 
one could consider a method of correcting the formal solution at 
these points, but the corrections must then also be applicable to 
the explicit construction of the A-operator 

A positive generalised opacity is assured by the following 
definition: 



Xi=Xi+^ aip\. 



(7) 



which only incorporates the physical opacity and a purely pos- 
itive contribution from the velocity field. The remaining term 
4fl//; must then be incorporated into the formal solution of the 
equation of radiative transfer. This term then has the same effect 
as if it was included in the opacity definition. For positive a/, the 
contribution to the formal solution is negative and therefore de- 
creases the intensity along the ray in a way similar to an increase 
in opacity provides. For negative a/, the contribution is positive 
which behaves like a negative opacity. This new solution indeed 
provides identical results to the previous discretisation method 
(see Sect.[3]l. 

After integrating Eq. |4]along a photon path from path length 
s\ to S2, the formal solution between two spatial points on a char- 
acteristic can then be written in terms of optical depth, since the 
optical depth t,\ along the characteristic relates to the path length 
by means of dr/ = ;f';di. 



(8) 



(9) 



//(t2) = //(TiK'-^^+ ri,e^-^MT- rSie^-'^dT 

where r,- = Ti(si) and 

Si = ^(s,-^-(pJI,^i+pIIm)] 
XI \ XI ^ 7 

Si ^ ^ ([1 - ^] [pjli-i + pIIm) + [4 + [1 - ^] p[] /,) (10) 
Xi 

In a discrete atmosphere model, the integrals in Eq.[8]can be 
evaluated by means of piecewise linear- or parabolic interpola- 
tion of the auxiliary source functions: 



r 
r 



Sie^ '"'dr ^ OijS i-ij + Pi^iS ij + yijS i+x^i 



Sie^ '"'dr ^ aijS i-ij + fSijS ij 



(11) 
(12) 



See lChen et al ]( l2007h for a more complete description. 



where the coe f ficient s a, 0, and y are described in 
lOlson & Kunasj d 19871) and iHauschildtl (|1992|) . Note that the 
contribution from the explicit wavelength derivative has been 
only linearly interpolated, whereas the implicit part can also be 
parabolically interpolated. 
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(a) 



(b) 



Fig. 1: Comparison of the comoving spectra for the cases ^ = 
in (a) and ^ = 1 in (b) where the new solution is represented by 
the black dashed line and the old one by the solid grey line. The 
wavelength scale is given in km/s around a fictitious line center 
The line is omitted and the continuum is purely absorptive. 



The construction of the matrix equation for the formal solu- 
tion as well as the analytic construc t ion sc heme for a A opera- 
tor described in lBaron & Hauschildtl (|2004| ) can be used without 
modification in the new formal solution. These steps will there- 
fore not be repeated here. 

3. The test calculations 

To compare the new solution with other solutions , we us ed a 
simplified test set up (see also iHauschildt & BaronI (|2004 and 
iKnop et al] (|2007|) ). The opacity consists of a grey continuum 
and a single atomic line. The continuum opacity Xk is grey and 
varies with the radial structure according to Xk °^ r^^- The ex- 
pression = K^ + cTii includes absorption as well as scattering 
(Tk, which are related by a thermalisation parameter - k^Ixk- 
The corresponding opacities /fune, criine, and thermalisation eune 
are used for the — assumed to be Gaussian shaped — spectral 
line of the two-level-atom, where its strength is parameterised 
relative to the continuum opacity R - xime/XK- The atmosphere 
has an extension of r,ni„ < r < 101 , where rmin - lO'^cm. A 
radial optical-depth scale is mapped onto the radial grid with the 
continuum opacity and ranges from 10"* < r < 10"*. To include 
wavelength couplings in the atmosphere, a velocity field is im- 
posed on the radial structure. 

At first, we checked that, for an absent velocity field, the old 
and new formal solutions provide identical results. Then we used 
monotonic velocity fields, in order to compare the new formal 
solutio ns with the well tested solution described in iHauschildi 
( Il992h . For the velocity field, a linear law v(r) = Vmax '"max/'" 
with Vtnax = lOOOkm/s was assumed. 

In the upper panel of Fig.[Tal the comoving frame spectrum 
for a pure absorptive continuum without a line — 0-^ = and 
R - — is shown for ^ = 0. The new solution is plotted with 
a dashed black line, while the old solution is represented by a 
solid grey line. The lower panel of Fig. [Ta] illustrates the dif- 
ference between the two solutions, implying that both solutions 
reproduce the expected flat continuum and differ about 1%. For 
comparison, the computations in Fig. [Ta|were repeated for the 
case ^ - I and are shown accordingly in Fig. [lb] The fully 
implicit ^ - I wavelength discretisation had difficulties in re- 
producing a flat continuum because to its diffusive behaviour, as 
also demonstrated by .Hauschildt & BaronI (l2004l) . This is also 
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Fig. 2: Same as in Fig.lT] apart from the fact that the spectral line 
is omitted and the continuum has a thermalisation parameter of 
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(a) 



(b) 



Fig. 3: Same as in Fig. [T] apart from the fact that the spectral 
line is strongly scattering (eiine = 10"^) and the continuum has a 
thermalisation parameter of e,^ - 10"'. 



true for the new solution, as can be seen in the upper panel of 
Fig. [TbJ the new formal solution closely follows the old solu- 
tion, as is also evident in the lower panel of Fig. [lb] The relative 
difference never exceeds 1 %, which is less than the deviation of 
the continuum from flatness. 

Since no scattering was included in the opacity, there is no 
need for an ALI and a single A-step will solve directly the sys- 
tem. The comparisons between the comoving spectra then di- 
rectly yield the differences in the formal solutions. 

In the next test, a scattering continuum opacity — of 6,^ = 
10^2 — was added, but the line was still omitted. The results for 
the case of ^ = are shown in Fig. |2a] and for ^ = 1 in Fig. |2b] 
with the same colour codes as before. The top panels show the 
comoving spectra and the relative differences are shown in the 
bottom panels. The differences between the scattering continua 
are of the order of only 10"^ %. The diffusive behaviour of the 
fully implicit solutions appears again in the deviations from a 
flat continuum. 

In further tests, a strong scattering line — with 7? = 100 and 
fiine = 10""^ — was added to the continuum opacity which is also 
affected by additional scattering — = 10"'. The results for 
the ^ = discretisation are shown in Fig. |3a] The upper panel 
again shows the comoving spectrum with the new solution in 
dashed black and the old solution in solid grey. The resulting line 
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Fig. 4: Alternating velocity field in terms of optical depth. 

profile is a typical scattering line in an expanding atmosphere. 
The apparent excellent match of the two solutions is verified by 
the plot of the relative differences in the lower panel of Fig. [3al 
The maximal difference is of the order of 10"^ %. 

The same calculation has been performed for the ^ = 1 
weighted wavelength discretisation. The results are shown in 
Fig. [3b] with the same colour codes as before. The comoving 
spectra show a disturbed continuum in the transition region be- 
tween the line and the continuum, which can again be attributed 
to the more diffusive behaviour of this wavelength discretisation. 
However, again the match between the new and the old solution 
is excellent. 

The excellent agreement between the old and new comoving 
spectra in the presence of scattering — with relative differences 
of the order of 10"-' % — may at first seem unusual because 
the differences in the formal solution itself were of the order 
only of one percent (see Figs. [Ta| and [Tbl l. However, in these 
cases, no ALI was needed due to the lack of scattering. In the 
scattering cases, an ALI was performed. Since only exact ele- 
ments from the A-operator have been used in the construction 
of the A* -operator, the operators derived from the old and new 
formal solutions also differ. From the excellent agreement of the 
final results, it is obvious that these differences cancel each other 
and enable far closer agreement for the moments of the radiation 
field in the ALI than the formal solution for a given source func- 
tion. 

In the preceding tests, there was no need to employ the new 
solution because the velocity fields were monotonic. To test 
the new solution for alternating wavelength couplings, a non- 
monotonic, alternating velocity field was used. 

To demonstrate the robustness of the new solution, we per- 
formed calculations for an alternating velocity profile — shown 
in Fig. |4] — that penetrates into the deep layers of the atmo- 
sphere. This run of the velocity profile was not even remotely 
possible within the old framework, and hence no comparison 
with the old solution is shown. The comoving spectrum is shown 
in Fig. |5] The peculiar features and broadness of the line result 
from the extreme velocity fields in the atmosphere. The ^ = 1 
discretisation had to be used, since the formal solution was ill- 
conditioned in the ^ = case. 

This calculation was not intended to model a physical atmo- 
sphere. It instead was intended to demonstrate the robustness of 
the new solution for even extreme wavelength-coupling condi- 
tions. 

We note that the removal of the 4a/// term from the opacity 
and its addition as an explicit source term has apparently little 
influence on the solution. In all of our test calculations, all dif- 
ferences with respect to established solutions were minor In the 
case of scattering in the line as well as in the continuum the dif- 
ferences become almost negligible. In comparison, the different 




-40000 -30000 20000 40000 

distance from linecenter [km/s] 

Fig. 5: Co-moving spectrum for the presence of the non- 
monotonic velocity field shown in Fig. |4] for ^ - I- The wave- 
length scale is given in km/s around a fictitious line center. The 
line is strongly scattering (eune = 10"'*) and the continuum has a 
thermalisation parameter of - 10"^. 

discretisations had a far more dominant effect on the solution 
than the addition of the explicit term, as can be seen in the case 
of the ^ = 1 weighting for which the diffusive behaviour is traced 
out. We should note further that the speed of convergence of the 
formal solution did not change significantly. 

The new solution method appears to be well suited to general 
radiative-transfer modelling involving non-monotonic velocity 
fields. This has immediate applications to ID spherical symmet- 
ric atmosphere models, e.g. for pulsating atmospheres of vari- 
able stars as well as shock fronts in supernova flows. 

It remains to be seen whether the different numerical nature 
of the solution, due to an additional explicit term, will be trou- 
blesome in some cases. Our test computations did not identify 
an example of numerically unstable behaviour during the explo- 
ration of the parameter space. 

4. Conclusions 

We have presented a new formal solution of the equation of ra- 
diative transfer. Its advantage with respect to established solu- 
tions lies in the avoidance of negative opacities that occur by de- 
sign in non-monotonic velocity fields. It allows the computation 
of the arbitrarily wavelength-coupled, radiative-transfer equa- 
tion that arises with complete avoidance of negative, generalised 
opacities that render the solution unstable. The new approach is 
required for any modelling that involves alternating wavelength 
couplings. Since one considers radiative transfer in multidimen- 
sional situations, non-monotonic velocity fields do naturally oc- 
cur. 

The relative differences between the new and old solutions 
was not found to exceed 1 % in all test calculations and were 
in most cases several orders of magnitude smaller Therefore, 
the new solution is a viable replacement of the old solutions and 
useful in modelling the atmospheres of variable stars or radiative 
transfer in hydrodynamically obtained structures. 
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